#include "field.hpp"
#include "meshTypeFunctions.hpp"
#include "bc_constantWaterDepth_velocity.hpp"
#include "util.hpp"

#include <catch2/catch_test_macros.hpp>
#include <catch2/catch_approx.hpp>

using namespace hms;
using Type = hms::boundaryCondition::Type;
using Catch::Approx;

TEST_CASE( "ConstantWaterDepth_velocity" ){
	auto test = []( auto&& dummyMesh ){
		using MeshType = remove_qualifiers<decltype(dummyMesh)>;
		auto mesh { uniMeshAs<MeshType>(4, 4, 1, 1, true) };
		
		using BC = ConstantWaterDepth_velocity<MeshType>;

		const BoundaryPatch& left   { mesh.boundaryPatch("left") };
		const BoundaryPatch& right  { mesh.boundaryPatch("right") };
		const BoundaryPatch& bottom { mesh.boundaryPatch("bottom") };
		const BoundaryPatch& top    { mesh.boundaryPatch("top") };

		Field depth { mesh, 1 };
		Field velocity { mesh, 2 };

		auto setBC = [&]( const BoundaryPatch& p, scalar dGhostArg ){
			velocity
				.template setBoundaryCondition<BC>( p, depth, dGhostArg )
				.update();
		};

		auto setBCs = [&]( scalar dGhostArg ){
			for ( const auto& p : mesh.boundaryPatches() ){
				setBC( p, dGhostArg );
			}
		};

		// subcritical outflow
		depth.values() = 1.0;
		velocity.values() = 1.0;

		scalar
			dGhost {1.0},
			dInner {1.5},
			vInnerX {1.0},
			vInnerY {1.0};

		depth.domainValues(left) = dInner;
		depth.domainValues(right) = dInner;
		depth.domainValues(top) = dInner;
		depth.domainValues(bottom) = dInner;

		setBCs(dGhost);

		auto wavespeed = [](scalar dInner, scalar dGhost){
			return 2*( sqrt(phys::g*dInner) - sqrt(phys::g*dGhost) );
		};

		REQUIRE(
			velocity.boundaryCondition(left).type()
				== Type::constantWaterDepth_velocity
		);
		
		/* dInner > dGhost, therefore the wave speeds are always in outwards
		 * direction */
		REQUIRE ( (velocity.ghostValues(left).row(0) == vInnerX - wavespeed(dInner, dGhost) ).all() );
		REQUIRE ( (velocity.ghostValues(left).row(1) == vInnerY).all() );

		REQUIRE ( (velocity.ghostValues(right).row(0) == vInnerX + wavespeed(dInner, dGhost) ).all() );
		REQUIRE ( (velocity.ghostValues(right).row(1) == vInnerY).all() );

		REQUIRE ( (velocity.ghostValues(top).row(0) == vInnerX).all() );
		REQUIRE ( (velocity.ghostValues(top).row(1) == vInnerY + wavespeed(dInner, dGhost) ).all() );

		REQUIRE ( (velocity.ghostValues(bottom).row(0) == vInnerX).all() );
		REQUIRE ( (velocity.ghostValues(bottom).row(1) == vInnerY - wavespeed(dInner, dGhost) ).all() );

		// supercritical outflow over boundaries
		depth.values() = 1.0;
		velocity.values() = 1.0;
		dGhost = 1.0;
		vInnerX = 3.5;
		vInnerY = 1.0;

		velocity.domainValues(left  ).row(0) = -vInnerX;
		velocity.domainValues(right ).row(0) =  vInnerX;
		velocity.domainValues(bottom).row(1) = -vInnerX;
		velocity.domainValues(top   ).row(1) =  vInnerX;

		setBCs(dGhost);
		
		REQUIRE( (velocity.ghostValues(left)   == velocity.domainValues(left)).all() );
		REQUIRE( (velocity.ghostValues(right)  == velocity.domainValues(right)).all() );
		REQUIRE( (velocity.ghostValues(top)    == velocity.domainValues(top)).all() );
		REQUIRE( (velocity.ghostValues(bottom) == velocity.domainValues(bottom)).all() );

		// sub- & supercritical flow over boundaries
		scalar vLarge {3.5}, vSmall {1};

		dGhost = 0.5;
		dInner = 1.0;

		depth.values() = dInner;
		velocity.values() = vSmall;

		depth.ghostValues(left  ) = dGhost;
		depth.ghostValues(right ) = dGhost;
		depth.ghostValues(top   ) = dGhost;
		depth.ghostValues(bottom) = dGhost;

		velocity.domainValues(left  ).row(0) = -vSmall;
		velocity.domainValues(left  ).row(0).middleCols(1,2) = -vLarge;
		setBC( left, dGhost );

		REQUIRE( (velocity.ghostValues(left).row(0).middleCols(1,2)
			== velocity.domainValues(left).row(0).middleCols(1,2) ).all()
		);

		// left corner cells
		REQUIRE( velocity.ghostValues(left).row(0)(0) == -vSmall - wavespeed(dInner, dGhost) );
		REQUIRE( velocity.ghostValues(left).row(1)(0) ==  vSmall );
		REQUIRE( velocity.ghostValues(left).row(0)(3) == -vSmall - wavespeed(dInner, dGhost) );
		REQUIRE( velocity.ghostValues(left).row(1)(3) ==  vSmall );

		depth.values() = dInner;
		velocity.values() = vSmall;
		velocity.domainValues(right ).row(0) =  vSmall;
		velocity.domainValues(right ).row(0).middleCols(1,2) =  vLarge;
		setBC( right, dGhost );

		REQUIRE( (velocity.ghostValues(right).row(0).middleCols(1,2)
			== velocity.domainValues(right).row(0).middleCols(1,2) ).all()
		);

		// right corner cells
		REQUIRE( velocity.ghostValues(right).row(0)(0) == vSmall + wavespeed(dInner, dGhost) );
		REQUIRE( velocity.ghostValues(right).row(1)(0) == vSmall );
		REQUIRE( velocity.ghostValues(right).row(0)(3) == vSmall + wavespeed(dInner, dGhost) );
		REQUIRE( velocity.ghostValues(right).row(1)(3) == vSmall );

		depth.values() = dInner;
		velocity.values() = vSmall;
		velocity.domainValues(top   ).row(1) =  vSmall;
		velocity.domainValues(top   ).row(1).middleCols(1,2) =  vLarge;
		setBC( top, dGhost );

		REQUIRE( (velocity.ghostValues(top).row(1).middleCols(1,2)
			== velocity.domainValues(top).row(1).middleCols(1,2) ).all()
		);

		// top corner cells
		REQUIRE( velocity.ghostValues(top).row(1)(0) == vSmall + wavespeed(dInner, dGhost) );
		REQUIRE( velocity.ghostValues(top).row(0)(0) == vSmall );
		REQUIRE( velocity.ghostValues(top).row(1)(3) == vSmall + wavespeed(dInner, dGhost) );
		REQUIRE( velocity.ghostValues(top).row(0)(3) == vSmall );

		depth.values() = dInner;
		velocity.values() = vSmall;
		velocity.domainValues(bottom).row(1) = -vSmall;
		velocity.domainValues(bottom).row(1).middleCols(1,2) = -vLarge;
		setBC( bottom, dGhost );

		REQUIRE( (velocity.ghostValues(bottom).row(0).middleCols(1,2)
			== velocity.domainValues(bottom).row(0).middleCols(1,2) ).all() 
		);

		// bottom corner cells
		REQUIRE( velocity.ghostValues(bottom).row(1)(0) == -vSmall - wavespeed(dInner, dGhost) );
		REQUIRE( velocity.ghostValues(bottom).row(0)(0) == vSmall );
		REQUIRE( velocity.ghostValues(bottom).row(1)(3) == -vSmall - wavespeed(dInner, dGhost) );
		REQUIRE( velocity.ghostValues(bottom).row(0)(3) == vSmall );
	};

	SECTION("UniMesh"){
		test(UniMesh{});
	}
	SECTION("RectMesh"){
		test(RectMesh{});
	}
	SECTION("StructMesh"){
		test(StructMesh{});
	}
	SECTION("Mesh"){
		test(Mesh{});
	}
}